function modelEq = projectionStepEquityPrice(model,setup,consClaim)

% Save output for the core of the model to compute equity price
thetaTrans = model.theta';
setup.out  = saveStatesControls(thetaTrans(:),setup);

% Adding loadings for dividends to out
if consClaim == 1
    posDiv        = find(strcmp(model.labely,'$c_t$'));   % defining equity on a consumption claim
else
    posDiv        = find(strcmp(model.labely,'$div_t$'));  % defining equity on a dividend claim from firms
end
setup.out.d0  = model.g0(posDiv,1);
setup.out.dx  = model.gx(posDiv,:);
setup.out.dxx = model.gxx(posDiv,:,:);
setup.out.dxxx= model.gxxx(posDiv,:,:,:);

% Settings for the optimizer
lambda = 0.01;
deltaOption = 1*0;

% Computing the equity price
% Starting values: use consumption
posC      = find(strcmp(model.labely,'$c_t$'));
nx = setup.nx;
if setup.appOrder == 1  
    thetaEquityPrice = zeros(1,1+nx);
elseif setup.appOrder == 2
    thetaEquityPrice = zeros(1,1+nx+nx*(nx+1)/2);
elseif setup.appOrder == 3
    thetaEquityPrice = zeros(1,1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6);
end
thetaEquityPrice0(1,:) = model.theta(posC,:);

setup.dispOn = 0;   %no display of optimization steps.

setup.nDim1Theta = size(thetaEquityPrice,1);
setup.nDim2Theta = size(thetaEquityPrice,2);
setup.scaleX_in_g = setup.scaleX_in_g*0+1;  %no scaling as expression for div is NOT scaled
if consClaim == 1
    if isfield(model,'thetaEquityPriceCons')
        thetaEquityPrice = model.thetaEquityPriceCons;
    else
        thetaEquityPrice = Projection_LMoptimizer(thetaEquityPrice0(:),setup,...
            setup.tolF,setup.MaxIter,lambda,deltaOption,setup.dispOn);
    end
else
    if isfield(model,'thetaEquityPrice')
        thetaEquityPrice = model.thetaEquityPrice;
    else
        thetaEquityPrice = Projection_LMoptimizer(thetaEquityPrice0(:),setup,...
            setup.tolF,setup.MaxIter,lambda,deltaOption,setup.dispOn);
    end
end
[~,resEuler] = NLSobjectFunc(thetaEquityPrice,setup);
MAEequity   = mean(abs(resEuler(:)));
RMSEequity  = sqrt(mean(resEuler(:).^2));
MaxAEequity = max(abs(resEuler(:)));

% Saving equity loadings
params         = model.params;
DIVss          = exp(model.gSS(posDiv,1));
tmp            = params.MUZss^(1-params.CHI*(1-params.CHI0)-params.CHI0);
EQss           = params.BETTA*tmp*DIVss/(1-params.BETTA*tmp);
eqSS           = log(EQss);
eq0            = thetaEquityPrice(1,1);
eqx            = thetaEquityPrice(1+1:1+nx,1)';
eqxx           = zeros(1,nx,nx);
eqxxx          = zeros(1,nx,nx,nx);
if setup.appOrder > 1
    idx   = 0;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            eqxx(:,i,j)    = thetaEquityPrice(1+nx+idx,1);
            eqxx(:,j,i)    = eqxx(:,i,j);
        end
    end
end
if setup.appOrder > 2
    idx   = 0;
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                idx = idx + 1;
                eqxxx(:,alfa1,alfa2,alfa3)    = thetaEquityPrice(1+nx+nx*(nx+1)/2+idx,1);
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %eqxxx(:,alfa1,alfa1,alfa3)= eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa1,alfa3,alfa1) = eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa3,alfa1,alfa1) = eqxxx(:,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %eqxxx(:,alfa1,alfa2,alfa2)= eqxxx(:,alfa1,alfa2,alfa3);                
                    eqxxx(:,alfa2,alfa1,alfa2) = eqxxx(:,alfa1,alfa2,alfa3);                
                    eqxxx(:,alfa2,alfa2,alfa1) = eqxxx(:,alfa1,alfa2,alfa3);                                
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %eqxxx(:,alfa1,alfa2,alfa3) = eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa1,alfa3,alfa2) = eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa3,alfa1,alfa2) = eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa3,alfa2,alfa1) = eqxxx(:,alfa1,alfa2,alfa3);
                    eqxxx(:,alfa2,alfa3,alfa1) = eqxxx(:,alfa1,alfa2,alfa3); 
                    eqxxx(:,alfa2,alfa1,alfa3) = eqxxx(:,alfa1,alfa2,alfa3);                 
                end
            end
        end
    end
end

modelEq              = model;
if consClaim == 1
    modelEq.labely      = [model.labely,'$pc^m_t$'];
else
    modelEq.labely      = [model.labely,'$p^m_t$'];
end
modelEq.gSS         = [modelEq.gSS;eqSS];
modelEq.g0          = [modelEq.g0;eq0];
modelEq.gx          = [modelEq.gx;eqx];
modelEq.gxx         = [modelEq.gxx;eqxx];
modelEq.gxxx        = [modelEq.gxxx;eqxxx];
modelEq.ny          = modelEq.ny+1;
modelEq.Gxx         = reshape(modelEq.gxx,modelEq.ny,nx^2);
modelEq.Gxxx        = reshape(modelEq.gxxx,modelEq.ny,nx^3);
if consClaim == 1
    modelEq.MAEequityCons   = MAEequity;
    modelEq.RMSEequityCons  = RMSEequity;
    modelEq.MaxAEequityCons = MaxAEequity;
    modelEq.thetaEquityPriceCons = thetaEquityPrice;
else
    modelEq.MAEequity   = MAEequity;
    modelEq.RMSEequity  = RMSEequity;
    modelEq.MaxAEequity = MaxAEequity;
    modelEq.thetaEquityPrice = thetaEquityPrice;
end

